Stroboscopic wavepacket description of non-equilibrium many-electron problems 
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We introduce the construction of a orthogonal wavepacket basis set, using the concept of stroboscopic time 
propagation, tailored to the efficient description of non-equilibrium extended electronic systems. Thanks to 
three desirable properties of this basis, significant insight is provided into non-equilibrium processes (both 
time-dependent and steady-state), and reliable physical estimates of various many-electron quantities such as 
density, current and spin polarization can be obtained. The use of this novel tool is demonstrated for time- 
dependent switching-on of the bias in quantum transport, and new results are obtained for current-induced spin 
accumulation at the edge of a 2D doped semiconductor caused by edge-induced spin-orbit interaction. 
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Wavepackets (WP) are a very useful concept when analyz- 
ing quantum mechanical scattering processes, since they com- 
bine local and wave-like aspects on an equal footing. Some 
of their more recent applications range from studies of the 
intrinsic spin Hall effect in semiconductors jl] 0], spin-flip 
dynamicsy], thermal averaging and its influence on interfer- 
ence patternsJfH] or transport of an electron through Luttinger 
liquidl!]. However, the use of traditional WPs in degener- 
ate fermionic systems raises difficulties since the exclusion 
principle restricts the available eigenstates that are superposed 
within a single WP. Several orthogonal wavepacketfl H 
and wavelet[9] approches were put forward in the past to 
accomodate the exclusion principle; however in contrast to 
our WPs these do not directly relate to typical many-electron 
states such as the electronic ground state or moderate pertur- 
bations from it at zero temperature. 

If we forego the time-dependent feature of WPs, the lat- 
ter problem is conveniently resolved with the introduction of 
Wannier functionsl Tol 11]: by occupying a finite number of 
them, we locally recover the exact eigenstates of a system of 
non-interacting electrons. 

In this work we combine the advantages of Wannier func- 
tions for extended systems with the time-dependent descrip- 
tion of WP propagation. This is achieved by generalizing the 
orthogonal WPs introduced by Martin and Landauer lll2ll for 
ideal ID leads. Our wave-packet basis set (WPB) has follow- 
ing three properties: (1) each basis function (WP) is localized 
in space, (2) occupying a subset of the WPB we recover the 
exact non-interacting many-electron ground state of a refer- 
ence Hamiltonian, (3) the WPB is generated by time propaga- 
tion through successive time-steps, T, of an initial set of WPs, 
according to a reference Hamiltonian. 

From the above properties it follows that we can view 
the whole basis set as a stroboscopic pictures of a continu- 
ous time-evolution of a suitably chosen family of initial WPs 
(Fig. [TJ. Since all WPs are orthonormal, each copy can be oc- 
cupied by precisely one electron and in time T each electron 
will move into its neighboring' WP. Similarly, if a single elec- 
tron is in a superposition of several WPs, in time T it will be 



in the same superposition but of the WPs obtained from the 
former by a single shift of the basis functions. This picture is 
valid as long as the reference Hamiltonian is time-independent 
in the region where the concerned WPs are localized. We will 
refer to this region as the bulk and to the rest, typically a much 
smaller region than the bulk, as the scatterer. Similarly, the 
bulk (scattering) WPs are those WPs that are generated with 
the bulk (bulk+scatterer) Hamiltonian. 

To obtain the time-dependent dynamics in the scatterer one 
needs to perform a full time-dependent simulation of the bulk 
WPs entering the scatterer. After certain time, the scattering 
WPs will return into the bulk where those WPs can once again 
be expanded into the bulk WPB and propagated as moves of 
duration t between the bulk WPs, i.e. analytically. Hence, 
the WPB offers a very simple interpretation of the processes 
as well as a framework to perform numerical time-dependent 
simulations. 

The consistency of the conditions (1) and (3) demands that 
the reference Hamiltonian posses translational symmetry in 
the direction of propagation. Its eigenstates in the Bloch 
form will be sufficient to create a basis such that each WP 
from the initial set will be spatially localized and their time- 
propagated WPs will slowly disperse with increasing time. 




— V 


w 


E 








k 



FIG. 1: Squared amplitudes of orthogonal stroboscopic wavepackets 
obtained by time propagation of the initial WP (in the centre) by 
a constant time-step T. The right- (blue/full) and left- (red/dashed) 
WPs belonging from the same energy band (inset) are shown. These, 
together with the WPs coming from the bands covering the rest of 
the spectrum, form a complete orthogonal basis set. 
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This property can be satisfied only if the reference Hamilto- 
nian is just that of the bulk. We may also construct the WPB 
for the combined system where the reference Hamiltonian is 
that of bulk+scatterer, but the scattering will typically result 
in a strongly delocalized WP (e.g. transmitted and reflected 
components). However, the scattering- WPs can be easily ex- 
panded into the bulk WPB, a fact of which we will make use 
later. 

Definition of the WPB and its format properties. To de- 
fine the basis set let us take an extended system specified by 
the reference Hamiltonian H with a continuous spectrum of 
eigenenergies £ € (eo,°°), H\e,a) — e\e,a) . To each eigen 
energy we will generally have a set of degenerate single- 
particle eigenstates \e,a) , a = 1,2, ...,N e , forming all to- 
gether a complete orthogonal whose normalization we choose 
such that (e',a'\e,a) = 8(e- e')8 a a i. 

From the above set we can generate an orthogonal and com- 
plete wave-packet basis set (WPB) by first choosing the initial 
set of wave-packets [29] 



|ra,0,a) 




de'U a ^(e')\e',a'), n = 0,1,2, 



(1) 

for an arbitrarily chosen division of the spectrum into energy 
bands {(£„ , e „+i)}™=o , « = 1,2,..., N e with bandwidths 



£„ . The division into energy bands must cover 



Afc n — fc n+l 

the full spectrum of ft but otherwise can be chosen so as to suit 
the physical situation as discussed later. The unitary, energy- 
dependent matrix U aa i{e) represents the second freedom of 
choice in the construction of the WPB. In this work we will 
use U aa i(e) = 8 aa i which is satisfactory for our present pur- 
poses, but in general it can be used either to adopt the bulk 
WPB to the scattering processes involed i.e. by shifting WPs 
in cetrain bands (U aa > (e) = 8 aa ie~ lkxX °), or to improve the lo- 
calization of the WPs, in analogy with Wannier functions !!!! . 
All the functions {\n,Q,a)} n are orthogonal by definition, 
since they are linear combinations of eigenstates from disjunct 
energy bands. 

The construction of the WPB is completed by forward and 
backward time propagation of the initial set 



\n,m, a) 



»\n,0,a) 



= ±1,±2, 



(2) 



by regular, band-dependent time steps T™ = 2n/Ae". It is 
easy to verify that this choice of time step guarantees or- 
thonormality of consecutive wave-packets within each band 



(n,m,a\n,m! , Of) = 8 m 



(3) 



Due to the orthogonality of the WPs we can uniquely 
expand any eigenstate of the reference Hamiltonian into 
the WPB with expansion coefficients (e,a\n,m,a) = 
(Ae„)^ 1 / 2 exp{— ;emr„}, with £ G (£„,£„ + Ae„). Con- 
versely, combining this with EqsQ] and [2] one obtains that 
Y,m\ n , m , a ) {n,m,a\e,a) = |e,ot) , from which follows that 
the WPB is also complete since the original set of eigenstates 
is a complete one. 
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FIG. 2: Abrupt switching on of the bias in a ID wire. In response 
to the bias AV, the right-going WPs for x > (white, the previously 
unoccupied band) start to fill the WPs from the left and the occupied 
left-going WPs for* < become empty. The finite extent of each WP 
causes oscillations, with period T, of the resulting occupancy N m (t) 
and hecne the current measured at fixed x m = vpmt (below). 



It has been already pointed out that the division into bands 
can be exploited to optimize the basis set to the particular 
physical problem. A typical choice of the energy bands is 
to take £,f = Ep for a cetrain n and all a, where Ep is the 
Fermi energy of the system. This way the ground-state is de- 
scribed by occupying all of the WPs in the bands below Ep. 
This means that we need to consider only few WP or electrons 
even though we are describing the local ground state proper- 
ties of the infinite many-electron system exactly (see 12811 ). 
Similarly, the non-equilibrium state is obtained by imposing 
different effective Fermi energies for WPs with different val- 
ues of a. 

We will now demonstrate the use of the WPB on several ex- 
amples from two rapidly developing areas of condensed mat- 
ter physics - time-dependent and/or ab initio simulations in 
quantum transport, and spin accumulation due to spin-orbit 
coupling in 2D systems. 

Time-dependent quantum transport. Understanding the 
quantum transport of charge through nanoj unctions made of 
individual atoms or molecules will be essential for progress in 
nanoelectronics. Due to the short spatial scale and short times 
involved it is clear that transient phenomena play an important 
role in understanding the functionality of nanodevices. At the 
same time, it has been recognized that the correct treatment 
of interactions demandes a time-dependent formulation of the 
density- or current-density functional the ory! 13l |. While sev- 
eral exact methods have been put forward lll4ll5 , 3, . 



U, due 



to their inherent complexity, they give restricted insight into 
the processes involved. Here we show that the WPB can pro- 
vide this insight in an elegant fashion, as well as quantitative 
results for transient times, oscillations or steady-state current. 

As an example let us start with the a ID electron gas in 
which at time t — a finite potential difference is applied 
(Fig.|2]i. Anticipating the application of the bias AV, we split 
the occupied part of the spectrum of a Hamiltonian for free 
electrons into occupied bands to AV, AV to Ep — AV and 
Ep — AV to Ep (we will refer to the last band as n = a, i.e. 
the active band), and two unoccupied bands Ep to Ep + AV 
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and Ef + AV to °°. The energy-normalized eigenstates are the 
plane-waves (x\e, a) — e ,akx / \Jlnk, k — y/2e, and a = ± 
for right- and left- going states respectively. The resulting 
WPs (x|«,m,±), obtained according to the Eqs. 11121 are ex- 
amples of the bulk WPs mentioned above and are identical 
to the WPs employed by Martin and Landauer in thier anal- 
ysis of quantum noise lll2ll . Due to the hopping of electrons 
between the WPs in time t, the current at the position of the 
m-th WP carried by electrons in the active band is in general 
given as I(t) =N m (t) / T, where N m (t) is the occupation of the 
m-th WP. 

Switching on the bias AV at x = and t = will ener- 
getically align WPs from the highest occupied, band local- 
ized in x < 0, with the WPs from the lowest unoccpied band 
and localized in x > 0. A transient phenomenon for time 
t ~ 2n/Ef < x = 2n/AV, which needs to be analyzed by per- 
forming a time-dependent simulation, will be related to dy- 
namics of those occupied WPs that had for t < nonzero am- 
plitude for both x < and x > 0. After that the time-dependent 
many-electron dynamics for x > will result in a train of 
right-going scattering orthogonal WPs within the active band 



{x\a, l,+;t) = 



d£ 



Ef-AV \/2KkkV 



t ^ k y{kx-\^{t+H)) ^ (4) 



occupied for / = 0,-1,-2,..., where t(k) is the transmis- 
sion amplitude for the applied step potential. The oc- 
cupation N m (t) of the m-th bulk WP due to this train is 
N m (t) = 2£^ | (a,l,+;t |a,m,+) | 2 which, after substituting 
the above expressions, gives finally 

Im(t) = -[[ F dede't*(e')t(e)F t - mT (e'-e), (5) 

1 J JE F -AV 

where F t (a>) = (AV)~ 2 £ / t^exp{-/a)(/T- f)}- in fact, this 
result is equally valid for abrupt swithing in ID wire with an 
arbitrary scattering potential and it represents a generalisa- 
tion of the Landauer formula for non-linear time-dependent 
response to abrupt switching on. In the long time limit we 
have F t (co) — > 5((»)/(AV) and we recover the non-linear Lan- 
dauer formula / = J AV \ t(e)\ 2 de/K . 

More specifically for the ID wire case, we can put 
t (k) = 1 and perform the integration with the result N m {t) = 
4sin 2 (AVf/2)/(AV) 2 £ / -~{(Z-m)T + r}- 2 . In the Fig. we 
show this result, calculated by taking the first 10 terms of this 
series, i.e. accounting for 10 WPs, for which we get a well 
converged answer. The relaxation to the steady-state current 
is characterised by oscillations with period T, in agreement 
with calculations based on non-equilibrium Greens functions 
within a wide band model[15]. 

The WPB-based picture offers a natural framework for the 
memory-loss theorem 1 15, 1^1 stating the independence of the 
steady state on the transient changes in external potential. In- 
deed, from the moment when the potential attains its long- 
time static form, it takes only a finite time until the WPs ex- 
periencing the transient potential leave the scatterer into the 
bulk, never to return. After that the occupancies of all the 



WPs inside this region are determined by the scattering of the 
bulk WPs within the long-time static potential. 

Our treatment here also indicates that the WPB repre- 
sentation can be used to perform numerical ab initio time- 
dependent simulations within the TDDFT framework, i.e. ac- 
counting for time-dependent self-consistent field in the scat- 
tering region. The time-evolution of the bulk WP as they enter 
the scattering region needs to be done numerically, but as soon 
as the scattered WP leaves this region, by expanding it into 
few bulk WPs one can perform its time evolution algebraically 
in a closed form. The density, current density or any other 
many-electron property is obtained by summing contributions 
from all stroboscopic images of the propagated WP. While the 
WPs will typically extend over several atomic distances, rela- 
tively few of them will be needed to compute local properties 
close to the scattering reagion, i.e. for a jellium model of a 
sodium mono-atomic wire with one atom missing (creating a 
gap and hence depletion of charge and corresponding Friedel 
oscillations) is well converged to the exact density of an infi- 
nite system with the gap using about 20 occupied WPs 12811 . 
Detailed implementation of at the self-consistent mean-field 
(TD DFT) methodology will be reported elsewhere lfl8ll . 

Edge-induced spin Hall effect. It has been recently shown 
that the interplay between nonzero Rashba-Bytchkov spin- 
orbit (SO) coupling, the scattering off the edge and nonzero 
electric current along this edge leads to a universal spin polar- 
ization localized close to the edge of the 2D gas in GaAs quan- 



tum wells ll20l 12111 . In parallel, several other authors 11221 \23 



24] considered the spin-orbit (SO) coupling due to nonzero 
gradient in potential in-plane, Vso — ~ a E [& x ^^( r )] ' P\ 
where Cfe is the strength of the SO coupling, a is the oper- 
ator of spin, V(r) is the confining potential at the edge and p 
the momentum operator ll30ll . The edge-SO scattering, analo- 
gous to the mechanism behind impurity scattering in the bulk 
of the 2D gas, seems to lead to effects similar to the Rashba- 
Bytchkov mechanism. 

Both of these effects can be understood and analyzed within 
the WPB description, but here we concentrate on the edge- 
SO scattering. We consider a 2D electron gas confined in 
the xy(x > 0) half-plane, with its edge being described by a 
model potential V(r) = W9(— x) where 9 is the step func- 
tion. This model is appropriate for typical doping densities 
n ~ 10 12 cm~ 2 where the Fermi wavelength Xp ~ 20nm is 
much larger than atomic spacing, principally determining the 
abruptness of the edge. The current is imposed in the y di- 
rection. Fourier transforming y — > k y , the SO term takes the 
form Vso = (XeO z W 8(x)k y , i.e. electrons with up and down 
spins in the z direction experience different scattering poten- 
tial at the edge. For each k y we construct a WP, localized in 
the x direction and constructed from the eigenstates of a bulk 
2D electron gas. If we time-propagate an initial WPs with an 
average k x pointing towards the edge and identical for both up 
and down spin states (left-going WP), the reflected WPs for up 
and down spins will have two different phase shifts 0| n, and 
hence a mutual spatial shift Is with respect to one other. For 
the model described here the shift, calculated from scattering- 
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states' phase shift is 

ls = (-^(<!>i-<i>i)) = -4a E -8(2W-(e))al + 0(a 4 E ), (6) 

where the averaging is over the energy band of the considered 
WP and e = Q%+l£)/2. We know that WPs separated by the 
time-step T are orthogonal and we may place one electron in 
each WP. The non-equilibrium situation can be set in the stan- 
dard fashion: occupying the WPs with k y > up to Ef + AV 
and those WPs with k y < only up to Ef. Deep inside the 2D 
bulk this WPs' shift will not contribute to any spin polariza- 
tion because a series of occupied WPs within each band gives 
homogeneous density. However, since the up- and down-spin 
WPs are shifted, this shift must be directly related to the spin 
accumulation close to the edge so that to first order in ole 

f dk v , , 2oCf / 

"T - "I ~ / ^hn{k y ) ~ --±y/2E^AV, (7) 

Jocc £7* <* 



where n(k y ) = J2Ep — k^/n is the number of initial WPs 
with momentum k y . The dependence on the magnitude of the 
confinement, W comes only in the 3th order, which follows 
from Eq.|6]and|7] 

^("T-"l) = -^M^V, (8) 

and hence the actual magnitude of the confinement poten- 
tial is rather unimportant. Both of the results, Eqs. [7] and [8] 
agree very well with more involved and exact Green's func- 
tion based treatments which will be reported elsewhere lE5tl . 
and demonstrate the usefulness of the WPB concept not only 
for qualitative but also for reliable quantitative estimates. 

It is interesting to compare the edge-SO scattering with the 
Rashba-Bytchov mechanism. The latter gives ll2lll m — n\ = 
-a|(2£>)~ 3 / 2 AV /(127T 2 ), where a R is the strength of the 
Rashba coupling; in the 2D GaAs systems it attains values lf26ll 
a R ~ 1.8 x KT 10 eV cm = 1.55 x l(T 2 a.u.*. On the other 
hand, the estimates for (Xe in GaAs quantum wells give [27] 
(Xe ~ 5.3A 2 = 5.53 x 10~ 4 a.u.*. The smallness of both ofe and 
Ur justifies the lowest order expansions used above. Finally, 
taking for the Fermi energy, Ef = 36meV= 3.01a.u.* corre- 
sponding to densities n ~ 10 12 cm~ 2 we find that the Rasba- 
mechanism is three orders of magnitude smaller than the edge 
spin-orbit scattering. In principle this might change at very 
low densities since the Rashba-mechanism increases while the 
edge SO scattering decreases with decreasing the Fermi en- 
ergy (or density) but for such low densities the behavior will 
be dominated by localization and interactions effects. 

In conclusion, our stroboscopic wavepacket basis permits 
both physical understanding and quantitative predictions to be 
obtained for a variety of non-equilibrium processes in which 
an extended system of electrons is subject to time-evolution 
while being coupled to bulk reservoirs. The stroboscopic 
construction permits the time-evolution of the system to be 
described straightforwardly, while the energy-localisation of 



the wavepackets within precise energy bands ensures that the 
Pauli principle is properly respected in coupling to the reser- 
voirs. 
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